Exactly solvable dynamics of the Eigen and the Crow-Kimura models 
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We introduce a new way to study molecular evolution within well-established Hamilton- Jacobi 
formalism, showing that for a broad class of fitness landscapes it is possible to derive dynamics 
analytically within the 1 /JV-accuracy, where N is genome length. For smooth and monotonic fitness 
function this approach gives two dynamical phases: smooth dynamics, and discontinuous dynamics. 
The latter phase arises naturally with no explicite singular fitness function, counterintuitively. The 
Hamilton- Jacobi method yields straightforward analytical results for the models that utilize fitness 
as a function of Hamming distance from a reference genome sequence. We also show the way in 
which this method gives dynamical phase structure for multi-peak fitness. 

PACS numbers: 87.10.+e, 87.15.Aa, 87.23.Kg, 02.50.-r 



I. INTRODUCTION 



Genome dynamics is an important problem in popula- 
tion genetics [1-3] and in molecular evolution [4-9]. Many 
authors investigated dynamics of evolution [10-13]. The 
Crow-Kimura and the Eigen models are very popular in 
evolution theory, describing quite well population genet- 
ics, the RNA virus evolution, and artificial evolution of 
molecules. The Crow-Kimura model describes evolution- 
ary process where mutation and selection are two parallel 
processes and describes mutations during the life time. 
The Eigen model describes the case where mutations oc- 
cur during the birth of new viruses (molecules) and is 
quite realistic for the RNA virus evolution. While exact 
solution is known for a simple case of single-peak fitness 
[14-16], there has been no success thus far in calculat- 
ing exact dynamics for a general fitness landscape. As 
in molecular evolution, there are numerous attempts to 
solve this problem at least approximately [10-13]. The 
fact is that evolution models are very subtle mathemati- 
cal objects and approximate solutions often give mislead- 
ing or inadequate results, especially in dynamics. Find- 
ing exact dynamics for these two models is well-known 
to be still an open issue. In this article we introduce 
Hamilton- Jacobi equations (HJE) as a mean to resolve 
it. These equations have been already applied in evo- 
lution theory to investigate population genetics of virus 
evolution with a finite population [l7|]. In Ref. [ij} HJE 
were applied and solved approximately for linear fitness. 
Also, HJE had been utilized in Refs. [18-19] to derive 
exact steady-state solutions for evolution models with a 
general fitness. In this work we show that it is possi- 
ble to obtain exact dynamical solutions of the Hamilton- 
Jacobi equations for the models where fitness is defined 
in terms of the Hamming distance from a reference (wild) 



sequence. The possibility of having analytical solutions 
that give the dynamics in a closed form is an important 
breakthrough in the theory of biological evolution. It al- 
lows one the investigation of a plethora of evolutionary 
pathways within one consistent formalism. By mapping 
evolution model to Hamiltonian mechanics and looking 
at the corresponding potential, it is possible to derive 
phase structure of the dynamics when exact dynamics 
are unavailable by other means. We show here the way 
to precisely calculate the movement of the maximum of 
the distribution for the population originally localized at 
a fixed distance from a reference sequence. This article 
is organized as follows. In Sec. II we review the known 
results for the Crow-Kimura model, analyze its dynam- 
ics via HJE when population is initially localized at some 
Hamming distance from a reference sequence, and inves- 
tigate the case when originally population is uniformly 
distributed across the sequence space. In Sec. Ill we 
solve the dynamics of the Eigen model. Our results are 
discussed in Sec. IV. 



II. THE CROW-KIMURA MODEL 



Main known results 



The 2 N genome configuration sequences are defined as 
chains of N spins s n , 1 < n < N, that can take on only 
two values s n = ±1- The reference configuration has 
all spins +1. The Hamming distance between a given 
configuration and the reference configuration is — 
s n )/2 = N(l — m)/2, where m is an overlap. This model 
describes the dynamics of probability distribution. We 
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denote configuration i by Si = 
the system is specified by 2 N 



The state of 



relative frequencies Pi, 1 < 



2 



< 2 N : 



where 



dPj 
dt 



— Aij Pj Pi ^2 Pj r j > 



Aij — $ij r j "i" m i 



(1) 



Here is the rate of mutation from configuration Sj 
to a new configuration Si, and is the fitness. Two 
configuration states have a Hamming distance = (N— 
J2k s f s j)/2: and m « = — JoN. When dij = 1 then my = 
7o and my = for d^ > 1 4]. For index i, the set of 
values 1 < i < 2 N is equivalent to the collection of N 

, we define the mean 



spins s fc . Identifying fo(si...s N ) 
fitness R: 



(2) 



The model defined here by Eq. ([I} had been introduced in 
Ref. [3] to describe the Drosphilla's evolution in a multi- 
allele model with simultaneously present mutation and 
selection processes. Because this model describes genet- 
ics of diploid evolution in infinite population the random 
drift is necessarily absent. The diploid evolution model 
of Ref. [3] is described by an equation in analogy with 
Eq.([T]) except that are linear functions of pi. In the 
model of [4] our Eq.|T]) describes an infinite population 
asexual evolution when there are either many alleles in 
one locus or many loci with two alleles in each. The selec- 
tion and mutation processes are decoupled in Eq.|T]), i.e., 
our model describes selection and mutation as parallel 
processes. This is different to a well-known model intro- 
duced by Eigen d, Q , where it is assumed that mutations 
originate as replication errors on the occasion of repro- 
duction events. Nowadays the Eigen's model is widely 
applied to describe the virus evolution. The model of 
Ref. [H as well as the Eigen's model [j| [§] have been 
suggested as molecular evolution models. Both these 
"connected mutation-selection" schemes of Refs. Q 
and "parallel" /'decoupled" scheme of Ref. [i| are simi- 
lar, giving similar pictures of evolution with only a slight 
difference in dynamics (e.g., see Fig. 1 in Ref. 15]). The 
difference between the connected multi-selection scheme 
and the parallel mutation-selection scheme of this work 
becomes transparent when both models are treated by 
a quantum Hamiltonian approach 0, [Til ]: the parallel 
scheme is described in terms of Hermitian Hamiltonian 
and the connected scheme is described in terms of non- 
Hermitian Hamiltonian. 

A value of R in steady state (dPi/dt = 0) is the main 
target of theoretical investigations. One can calculate 
R as maximal eigenvalue of a matrix Aij @,|. Con- 
nection between the Crow-Kimura model and quantum 
mechanics has been established in Ref. Q , where matrix 
—Aij has been identified with the quantum Hamiltonian 
H for N interacting quantum spins. One can calculate 
the maximal eigenvalue of the operator —H [1, HB] as 



R = lim 

/3— oo 



ln?Ycxp[-/?g] 

' 



(3) 



N 

# = ™$>£-i) 
fe=i 



/oW. 



'NJ 



(4) 



where cr|! and cr£ are Pauli matrices acting on the spin 
in the fcth position [16| . We are interested in symmetric- 
fitness case with /o(si-..sat) = Nf(J2k=i s k/N). For a 
symmetric fitness function and permutation-symmetric 
initial distributions all configurations at the Hamming 
distance I from the reference sequence (selected with 
s n = 1,1 < n < N) have one value of probability so 
the probability of selecting the entire class of configura- 
tions is cpi. For symmetric fitness the mean fitness is 
calculated as in Refs. 

R 



— = k = max U(x), 

N -1<X<1 ^ ' 



u(x) = f(x) - 1 + vT 



(5) 



The maximum point of Eq.([5]) occurs at x = x c . It fol- 
lows from Eq.© that x c can be interpreted as "bulk 
magnetization" in analogy with other models of statisti- 
cal mechanics [H 0, [2(| : 



NTrexp[-(3H] 



lim 



0- 



Despite the lack of direct biological meaning, we need 
to find x c to calculate the mean fitness. For symmetric 
fitness function and permutation-invariant original distri- 
bution there is a set of differential equations for (N + 1) 
relative probabilities pi , < / < N [5j : 

dpi _ 
dt 

■) / 

Pi[Nf{\ --)-N] + (N-l+ l)p,_ x + (I + l) Pl+1 . (6) 

Probability of finding all configurations at the Hamming 
distance I is pi/ J2k P k - Mapping of the system of nonlin- 
ear equations (1) onto the system of linear equations (6) 
was calculated in Refs. [20-21]). In Eq.(6) we omit p_i 
and pn+i for I — and I = N, and set 70 = 1. In biolog- 
ical applications a magnetization-like measure of surplus 
or surface magnetization can be defined as 



E,(l-2//JV)p, 



(7) 



The main goal of this work is to calculate the dynamic 
of x m from given initial distribution. Having the value 
of x c it is possible to calculate the value of x m in steady 
state by solving: 



f(x v 



k. 



(8) 



Various interpretations of bulk magnetization x c and sur- 
face magnetization x m were analyzed in Ref. 0, [20| . In 
next sections we solve the model for the dynamics and 
determine explicit role of x c for various sub-phases in 
dynamics. 
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B. HJE for Crow-Kimura model 

As in Ref. [l8|, at a discrete x = 1 — 21 /N we use 
the ansatz pi(t) = p(x,t) ~ exp[Nu(x, t)]. Equation (6) 
can be then written as Hamilton-Jacobi equation for u = 
\np(x,t)/N (in [l8| we gave an equation for individual 
probabilities in the sequence): 



U(x) 



^+H(u',x) = 0, 



where u' — du/da 



1 + x 2u' , 1 x „-2u' 



H(u\x) = f(x)-l + —e 2u + , } 



(9) 



(10) 



where the domain of a; is— 1 < ar < 1, and the initial 
distribution is u(x, 0) = uq(x). Equation (9) describes 
a class of probabilities and the equation describing one 
sequence of probabilities was given in Ref. [HI]. In the 
limit of t — > oo the asymptotic solution of Eq. (9) is 



u(x, t; k) = kt + Uk{x), 



(11) 



where Uk(x) can be calculated from Eq.(9) [la] and 
the mean fitness is Nk. Function U(x) in Eq.(5) has 
a simple physical interpretation as potential, i.e., the 
minimum of —H(u,v) with respect to v at a fixed x: 
U(x) = min v [—H(v,x)]. It is well-known from mechan- 
ics that motion is possible on an interval when energy of 
the system is larger than potential U(x) inside this inter- 
val. In maximum-principle approach the largest eigen- 
value is identified with the mean fitness k. Similarly, 
— A; is the maximal energy of the Hamiltonian H(y,x) in 
Eq. (10). A realistic hypothesis would be to assume that 
the asymptotic solution u(x, t; k) is stable against pertur- 
bations only if k is calculated according to Eq. (5). It is 
possible to obtain more results even without solving the 
dynamics exactly. We know from physics that motion in 
potential that has a single minimum is drastically differ- 
ent from motion in potential with two or more minima. 
Therefore, when in Fig. 1 function U (x) changes from 
that depicted by the continuous line to that presented 
by the dashed line but for potential well U(x) that has 
two maxima and two minima near x = we should an- 
ticipate phase transition. Here, we focus on the fitness 
f(x) = cx 2 /2 [4] (the solid curve in Fig. 1 corresponds to 
c = 2). It results from Eq. (5) that in this case U(x) has 
two extrema located on the interval [—1; +1]: the mini- 
mum at x = Q, and the maximum at x = x m . To solve 
Eq.© subject to these initial data we use a standard 
procedure [23l . l2~ij by allowing to reduce the correspond- 
ing partial differential equation to a system of ordinary 
differential equations. Namely, consider the following set 
of equations: 

x = H v (x, v) = -(l + x) e 2v + (1 - x) e- 2v , 

v = -H x (x, v) = f'{x) + (e 2v - e- 2 ")/2, (12) 
u — v H v (x, v) — H (x, v) — vx + q, 




Figure 1: Function U(x) = f(x) + \/l — x 2 — 1 for f(x) — x 2 
(solid curve) and for f(x) = 4exp(— 8 + 8x) (dashed curve). 
For the latter there are two extrema where U'(x) = 0: the 
maximum at 0.9995 (it is too high in is not shown in the 
graphics) and the minimum at 0.497. 



subject to the following initial conditions: x(0) = xq, 
v(0) = vq(xq), u(0) = uq(xo). Here, v = du/dx, v (x) = 
itg(a;), and q = du/dt. The corresponding solution of 
Eq. (|12p in the (x, t)-space is called the characteristic of 
Eq. ©. Further, Eqs. © and (12]) imply q = 0. Along 
the characteristic x = x(t) and variable q is constant, 
so q is selected to parameterize these curves. Using the 
equation q = f(x) - 1 + (1 + x)/2e 2v + (1 - x)/2e- 2v , we 
transform the first equation in Eq. (fT2"]) into 



±2^] 



q+l- f{x)Y +x 2 



(13) 



Having the solution of the characteristic system given 
by Eq.(12), we can derive the solution of the original 
Eq.© [24| by integrating the equation u = vx + q. For 
biology applications it is important to know motions of 
distribution maxima. For the purpose of finding these 
motions consider the following initial distribution 



uq(x) = —a(x — xq) 2 



(14) 



It is relatively easy to derive relaxation formulae for large 
values of parameter a. We can calculate them directly 
from Eq. (13), using equation q(x*,t*) = f(x*) for the 
maximum point location x* . The maximum of the distri- 
bution moves along the branch of Eq. (fT3"|) that preserves 
the sign of xq. By integrating Eq. (fT3"|) along the charac- 
teristic through the point (x* , t*) and assuming that x(t) 
does not change its sign, we are getting 



sgnx 



V(/(a:*) + l-/(0) 2 +£ 2 -l 



(15) 



If at some point x\ the characteristic x(t) changes its 
direction the point X\ can be determined from the con- 
dition 



[f(x*) + l-f(x 1 )} 2 +. 



1 = 0. 



(16) 
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Figure 2: The dynamics of the maximum point x(t) for the 
Crow-Kimura model (f(x) — x 2 ) for different initial values xo 
in the distribution (14). The continuous curves are analytic 
results of Eqs.(15) and (17). The symbols are the results 
of numerical solutions of the Crow-Kimura model given by 
Eq.(6), where N = 1000. 



In the latter case the integrals should be summed up over 
the intervals (xq,Xi) and {x* ,x\). This summation gives 



^ = sgnx ( . 



V(/(^) + i-/(0) 2 + e 2 -i 
I V(/(^) + i-/(0) 2 + e 2 -i 



(17) 



Let Ti be such that for t < T x Ea.(|T5]) holds, and for 
t > T x Ea.lfTT)) holds. At T x we have the condition 



Ti = 



sgn xo 



2 -i 



, (18) 



where X x is a root of [j{Xi) + 1 — /(xq)] -r ^ 



Xq - i = 0. 

For the quadratic fitness /(a;) = cx 2 /2 with c > a 
selective phase exists at c > 1. Then, x m = 1 — - and 
x c = Vl c~ 2 [4] . When i — > oo the maximum converges 
to .x = x m . To define the dynamics of the maximum at 
~x c < xq < x c we use Eqs. (fl~5]) and (fT7|) . where 



xi = sgnx 



V^x* 2 + 2(c - 1) - 2[(c - l) 2 - c 2 x* 2 ]V2 



. In the region where x c < |xo| < 1 we use Eq.(15). To 
find Xi in accordance with Eq. tfTB)) we use 



Xx = sgnx 



2[l-(l-x 2 )V2] 



(19) 



Figure 2 shows the evolution of the maximum for c = 2 
for xo = 0, 0.1, 0.3, 0.7, 0.95. These results demonstrate 
the excellent agreement of analytic solutions given by 
Eqs.(15) and (17) with the results of the numerical in- 
tegration of Eq.(6). Note, Fig. 2 shows that for xo < x m 
the maximum moves initially away from the wild config- 
uration and returns to its neighborhood in later times. 
The minimal x*(t) is just X\. If x*(t) describes the po- 



sition of maxima then v(x*(t),t) 
Eqs.dEJ) give 



dv(x'*(t),t) _ 
dt 



— and 



dx*(t) 
dt 



-2x*(t) 



u xx {x*{t),t)' 



x(0) = x , (20) 



where u xx {x,t) — dv/dx. The motion of the maximum 
of the distribution cither towards the wild sequence or in 
the opposite direction depends on the sign of f'(x*(t)) + 
2x*(t)u"(x*(t),t). 



C. The flat original distribution 

When any of 2 N configurations is uniformly populated 
then the initial condition for the entire probability class, 

having probability ( ( jv(i+x)/2)!)f^' y ields 



l + x 1 
Mo (X) = — m — 



(21) 



Solution (21) has a peak at x = 0. Let us calculate 
threshold-time T2 such that for t < T2 the population 
peak is in the class of x = 0. Assuming that at the 
moment t* the maximum is at point x*, we solve Eq.(13) 
for the characteristic with end-point (x*.t*) and, thus, 
take q — f(x*). The related characteristic curve starts at 
the point x(0) = x* , passes through the point (xi,t*/2) 
(xi is computed from Eq.(16)), turns, and finally reaches 
the point (x*,t*). Thus, Eq. fTfH) gives 



t* = sgnx 1 



X* 



1 



(22) 



Now we take the limit as x* — > and find the threshold 
time T2. When /(x) = cx 2 /2 and c > 1 this time is 



T 2 = c 0S - 1 ( v /l-l/c)/ % /c — T. 



III. THE EIGEN MODEL 



(23) 



As shown in Refs. [5-6], for 2 N probabilities Pi there 
is a set of equations 



dr 



2 N 2 N 

i=i j=i 



(24) 



Elements Qij of the mutation matrix give the proba- 
bilities that an offspring of configuration j belongs to 
configuration i. In this model mutations are quanti- 
fied by Qij = q N - d ^)(l - q) d ^) and 7 = N(l - q), 
where exp[— 7] = q N is the probability of having ex- 
act copy, Tj = /(l — 21/N) is the fitness, and I is the 
Hamming distance of the jth configuration from the ref- 
erence configuration. The Hamming distance between 
configurations i and j (that have spins spins s l n and 
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s l n , respectively) is d(i,j) = £) n (l ~ s n s n)- Consider- 
ing again the (N + 1) Hamming-class probabilities p\ for 
pi = exp[Nu(x,t)]andx = 1 - 21/N, Eq.(24) of Ref. [H 
has been mapped onto the following equation 



du 
dt 



f{x)e 



7[ch(2ti')+Ksh(2ti')-l] 



(25) 



where r = tiV. Asymptotic solutions u(x,t;k) = kt + 
itfe(x) (fc is a mean fitness |25|) in the limit of t — > oo are 
as follows 



k = max L/(x), U(a;) = f(x) cxp( 7 [-l + y/l - x 2 ]p6) 

—1<X<1 

where x r and obtained from 



U'(x c ) = 0, f(x m ) = f(x c ) exp(- 7 [l - Vl-^])- 

When x c < |a;o| < 1 then for initial distribution given 
by Eq.(14) with a >> 1 the position of the maximum 

(t*, £*) is 



sgnx 



. /(a; )y(ln/(|}+ 7 ) 2 - 7 2(l-e 2 ) 



(28) 



For all other cases the solution is 

«cm TV* / 

f* 



sgnxo / 
2 V 



+ 



,o /(.T*)y(ln^+ 7 ) 2 - 7 2(l-£ 2 ) 



)■ 



/(**) y(ln^+ 7 ) 2 - 7 2(l-e 2 ) 



where xi can be calculated from the condition 



In 



/(**) 



+ 7 -y(l-^) = . 



(29) 
(30) 



Finally, for relaxation from the flat distribution we get: 

Xl 

t* = sgnx* / ^ ^=,(31) 



»• /(^)^(ln^g+ 7 ) 2 - 7 2 (l-e 2 ) 
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Figure 3: Dynamics of maximum density points x*(t*) for the 
flat initial distribution given by Eq.(14): (a) Crow-Kimura 
model where (i) f(x) — Sx, (ii) f(x) — x 2 , (iii) f(x) = x 2 + 
0.2x 4 , (iv) f(x) = 4exp(a; - 1), and f(x) = 4exp(-8[l - a;]) 
(dashed line); (b) Eigen model where 7 = 2 and (i) f(x) = 
2(x + 1), (ii) f{x) = x 2 , and (iii) f(x) = exp(4x). Contin- 
uous curves are the analytical results. The symbols are the 
solutions of numerical integration. 
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Figure 4: The dynamics of the mean fitness R(t) for the 
Crow-Kimura model (f(x) = x) for different initial values 
xo = 0.5 in the distribution (14). The symbols are the results 
of numerical solutions of the Crow-Kimura model given by 
Eq.(6), where N = 1000. The upper line is an approximate 
result by diffusion method, the lower line is our exact result. 



equation. This method is qualitatively similar to semi- 
classical methods, well-known in quantum mechanics. 
Our approach has an accuracy of 1/N, where N is 
genome length. There is straightforward connection 
between our current method and methods that utilize 
statistical-physics analogies with Ising spins. Specifically, 
two different sub-phases that have been determined with 
our method describe two different relaxation regimes 



R(t) 
0.5 

0.4 

0.3 



IV. DISCUSSION 

We have considered discrete-error classes in contin- 
uum approximation, replacing the system of equations 
for molecular evolution by a single Hamilton-Jacobi equa- 
tion. Dynamics have been obtained by solving this 



0.2 0.4 0.6 0.8 1 

Figure 5: The dynamics of the mean fitness R(t) for the 
Crow-Kimura model (f(x) = x 2 ) for different initial values 
xo = 0.5 in the distribution (14). The symbols are the results 
of numerical solutions of the Crow-Kimura model given by 
Eq.(6), where N — 1000. The upper line is an approximate 
result by diffusion method, the lower line is our exact result. 
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c 


1.1 


1.2 


1.3 


1.4 


1.5 


1.6 


T 2 


2.397 


1.791 


1.466 


1.252 


1.098 


0.980 


ti 


3.998 


2.572 


1.953 


1.591 


1.351 


1.177 



Table I: Comparison of ti, the result of [20J for the threshold 
time period in case of initially flat distribution, with T2, our 
exact result by Eq.(23) for Crow-Kimura model with f(x) = 
cx 2 /2. 



(i.e., Eqs.(15) and (17) for Crow-Kimura's model; and, 
Eqs. (28-30) for Eigen's model). These two relaxation 
regimes correspond exactly to two different magnetiza- 
tion values as discussed in Refs. [E B 0| • Singularities x c 
in relaxation periods correspond to bulk magnetization. 
Initially, when the entire virus population is in one ge- 
netic configuration that is closer to the wild configuration 
than the sequences with the same value of x c , the maxi- 
mum in the population distribution moves to the steady 
state x m . This is in analogy with surface magnetization. 
On the other hand, when the initial configuration is far 
away from x c , the maximum of the distribution moves 
away from the wild configuration in the initial phase and 
moves towards later phase. The single minimum 

at x = of the evolution potential U(x) (i.e., Eq.(5) for 
Kimura's model, and Eq.(26) for Eigen's model) gives 
smooth dynamics (see Eqs. (15), (17), Fig. 2, and Eqs. (28- 
30)). Equations (22) and (31) give the evolution from 
the original flat distribution in the Crow-Kimura and the 
Eigen models, respectively. These results are presented 
in Fig. 3 for several choices of fitness function. Analytical 
dynamics of maximum-density points x*(t*) is in excel- 
lent agreement with numerical solutions for the original 
formulation of these models. The second phase of the 
dynamics with a jump in the position of x*(t) (seen as 
the dashed line in Fig. 3a) is related to the presence a po- 
tential well (indicated by the dashed line in Fig. 1). Pre- 
liminary numerical studies of similar problems indicate 
the existence of a similar phase with a jump that does 



not require a potential well but a steep potential. The 
evolution dynamics is a highly non-trivial phenomenon. 
As we demonstrated in this work, even for monotonic and 
smooth fitness landscapes it is possible to have discontin- 
uous dynamics in analogy with the punctuated evolution 
of Ref . [26] or the shock waves of Ref. [27| ) . Such discon- 
tinuous dynamics for smooth fitness function has been 
also found in Ref. [28| , where the dynamic of the evolu- 
tion model was investigated numerically for four-valued 
spins. In the current article we suggest the analytical 
method to investigate discontinuous evolution for a gen- 
eral fitness case. In Ref. (20| an analytic approximation 
that would be accurate for large c have been suggested 
for the dynamic of Crow-Kimura model. In Table 1 we 
compare our exact result for Ti obtained from Eq.(23) 
with the corresponding expression derived by the method 
of Ref. (20| (by setting A = 1 in Eqs. (4) and (65) of 
Ref. [13]). Our method gives the full distribution, while 
the method of Ref. [20( gives the position of the distri- 
bution maximum. In summary, we considered HJE to 
obtain exact dynamics and used Hamiltonian mechanics 
for qualitative analysis of evolution models. Our results 
are valid for any analytic fitness function. The diffusion 
method of Refs. [10-13] is valid only near the maximum of 
distribution or for the case of weak selection, and yields 
inaccurate results when applied for long relaxation peri- 
ods or for calculating mean fitness. This yields the error 
greater than 50% after t = 0.2 (see Figs. 4-5). The HJE 
approach is self-consistent, with no need to use genome 
length (which is in contrast to Refs. [10-13]), and gives 
the dynamic with the 1/N- accuracy. 
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